// MLE for predicted Wt. proportion in the diet of pred j
DATA_SECTION
  // read in data from predWt2.dat
  init_int LWnobs
  init_int Nnobs
  init_int wtdatnobs
  init_int maxdatnobs
  init_int npars

  init_vector Lp(1,LWnobs)
  init_vector Wp(1,LWnobs)
  init_vector prey_size(1,Nnobs)
  init_vector Np(1,Nnobs)
  init_vector propN(1,Nnobs)
  init_vector Wt(1,Nnobs)
  init_vector pred_size(1,wtdatnobs)
  init_vector wt(1,wtdatnobs)
  init_vector G(1,wtdatnobs)
 LOCAL_CALCS
   for (int l=1;l<=wtdatnobs;l++)
   {
     if(G(l)==0)
     {
       G(l)=0.0001;
     }
     if(G(l)==1)
     {
       G(l)=0.9999;
     }
   }
 END_CALCS
  init_vector WtTot(1,wtdatnobs)
  init_vector juse(1,maxdatnobs)
  init_vector z(1,maxdatnobs)
  init_vector pars(1,npars)
  !! cout << pars<<endl;  // tracks progress....
  
  int l // pre-allocate k
  int i // ditto
  int j //ditto
  
INITIALIZATION_SECTION
  logalphap 1.89
  betap 0.09
  sigmaw 3.2
  beta0g -2.6
  beta1g 1.86
  beta2g 0.24
  beta3g 0.0018
  beta4g -0.66
  logphig -0.9
  q 3.4
  logalphaz 1.9
  betaz .17
  logshapez -1.7

PARAMETER_SECTION


	number bb
	vector aa(1,maxdatnobs)
	number shapez // pre-allocate k
	number phig   // ditto
	init_number logshapez // pre-allocate k
	init_number logphig // ditto
	
	init_number logalphap
	init_number betap
	init_number sigmaw
	init_number beta0g
	init_number beta1g
	init_number beta2g
	init_number beta3g
	init_number beta4g
	init_number q
	init_number logalphaz
	init_number betaz
	
	vector logWphat(1,LWnobs) // used to find params logalphaz and betaz
	vector Wphat(1,Nnobs)  // used to predict lambda but should be length of Nnobs
	vector epsw(1,LWnobs) // usedto fit logWphat and log(Wp), should be length of LWnobs
	// vector f(1,LWnobs)	// not used
	// vector NLLw(1,LWnobs)	// not used 
	vector logzmean(1,maxdatnobs)	// used to find logalphaz, betaz, and shapez
	vector zmean(1,maxdatnobs)	// used to find """"
	vector NLLz(1,maxdatnobs)	// correct
	vector limhat(1,maxdatnobs)	// used to find the upper limit of prey sizes= should be length of pred.size
	vector lambda(1,wtdatnobs)	// used to estimate proportion of weight , should be length of pred.size
	vector logitphi(1,Nnobs)	// correct
	vector phi(1,Nnobs)	//correct
	vector mu(1,wtdatnobs)	// correct
	vector logitmu(1,wtdatnobs)	// correct
	vector a(1,wtdatnobs)	//correct
	vector b(1,wtdatnobs)	//correct
	//vector Guse(1,wtdatnobs) // not used	
	vector NLLg(1,wtdatnobs)	// correct
	objective_function_value obj_fun;

//PRELIMINARY_CALCS_SECTION
  
PROCEDURE_SECTION
  shapez = mfexp(logshapez);
  phig   = mfexp(logphig);
//  calc_LW();
//  cal_lim();
//  cal_lambda();
//  calc_wtprop();

//FUNCTION calc_LW(logalphap betap)
  // calculate the lw weight relatioships
  logWphat = logalphap+log(Lp)*betap;
  epsw     = log(Wp)-logWphat;
  obj_fun += -0.5*norm2(epsw); // double check that this is correct - need to be log normal errors....
  //f=(norm2(epsw)); 
 // obj_fun +=LWnobs/2.*log(f);    // make it a likelihood function?

//FUNCTION calc_lim(logalphaz betaz shapez)
  //calulate the max size prey available for each pred size - find logalphaz and beta z
  logzmean = logalphaz+betaz*log(juse);
  zmean    = mfexp(logzmean);
  bb       = shapez;  // scalar
  aa       = shapez*zmean; // vector
  //for (i=1;i<=maxdatnobs;i++)
  // neg log lik for gamma dist
  NLLz     = -1*(aa*log(bb)-gammln(aa)+elem_prod((aa-1),log(z))+-bb*z);
  obj_fun +=sum(NLLz);
  
//FUNCTION calc_lambda(logalphaz betaz logalphap betap q)
  for (j=1;j<=wtdatnobs;j++)
  {
    limhat(j) = mfexp(logalphaz+betaz*log(pred_size(j))); // max prey size predator of size j can consume
    // Wphat     = mfexp(logalphap)*pow(prey_size,betap);  // vector of predicted prey weights  should be length of Nnobs (not to be confused with logWphat above which is length of LWnobs)
   // logitphi  = elem_div(limhat(j)*q.-prey_size,prey_size); // vector of phis (logistic curve for each size bin of prey - should be length of Nnobs
   // phi       = elem_div(mfexp(logitphi),1+mfexp(logitphi)); // vector (length (Nnobs) convert with logit - phi is now on scale of 0-1.
	
	for (i=1;i<=Nnobs;i++) 
    	{
		phi(i)=1;
		low=limhat(j)-q
		if (prey_size(i)<low)
		{
			phi=1-(pow((prey_size(i)-low),2)/(low+pow((prey_size(i)-low),2)))
		}
		if (prey_size(i)>limhat(j))
		{
			phi=1-(pow((prey_size(i)-limhat(j)),2)/(limhat(j)+pow((prey_size(i)-limhat(j)),2)))
		}
//		if ( prey_size(i)==0)
//		{
//			phi(i)=1;
//		}
	}
	
    //lambda(j)=sum(elem_prod(elem_prod(Np,Wphat),phi))/sum(elem_prod(Np,Wphat));    // scalar based on biomass
    lambda(j) = sum(elem_prod(Np,phi))/sum(Np); // should be a scalar, Np = # of prey, at each size (prey_size)         # numerically based
  }
  
//FUNCTION calc_wtprop(beta0g beta1g beta2g beta3g beta4g phig)
  logitmu = beta0g+beta1g*lambda+beta2g*pred_size+beta3g*elem_prod(pred_size,pred_size)+beta4g*elem_prod(pred_size,lambda);
  mu      = elem_div(mfexp(mu),1+mfexp(mu));
  a       = mu*phig; // vector
  b       = (1.-mu)*phig;  // vector
  // neg log lik for beta dist
  NLLg=-1*(gammln(a+b)-gammln(a)+gammln(b)+elem_prod(a-1,log(G))+elem_prod(b-1,log(1.-G)));
  // NLL.g[is.na(NLL.g)]=1000000 penalize function?
  obj_fun +=sum(NLLg);
  
  //sum(sum.na(NLL.z),sum.na(NLL.w),sum.na(NLL.g))
 
REPORT_SECTION
  cout <<endl<<"========End of phase: "<<current_phase()<<" ============"<<endl<<endl;

  report <<"predicted proportion of wt in the stomach"<<endl;
  report <<mu<<endl;

TOP_OF_MAIN_SECTION
  arrmblsize = 10000000; 
